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Abstract 



We introduce a new model for describing the fluctuations of a tick-by-tick single 
asset price. Our model is based on Markov renewal processes. We consider a point 
process associated to the timestamps of the price jumps, and marks associated to price 
increments. By modeling the marks with a suitable Markov chain, we can reproduce 
the strong mean-reversion of price returns known as microstructure noise. Moreover, 
by using Markov renewal processes, we can model the presence of spikes in intensity of 
market activity, i.e. the volatility clustering, and consider dependence between price 
increments and jump times. We also provide simple parametric and nonparametric 
statistical procedures for the estimation of our model. We obtain closed-form formula 
for the mean signature plot, and show the diffusive behavior of our model at large scale 
limit. We illustrate our results by numerical simulations, and find that our model is 
consistent with empirical data on the Euribor future, [j] 

Keywords: Microstructure noise, Markov renewal process, Signature plot, Scaling limit. 



1 Tick-by-tick observation, from 10:00 to 14:00 during 2010, on the front future contract Euribor3m 



1 Introduction 



The modeling of tick-by-tick asset price attracted a growing interest in the statistical and 
quantitative finance literature with the availability of high frequency data. It is basically 
split into two categories, according to the philosophy guiding the modeling: 

(i) The macro-to-microscopic (or econometric) approach, see e.g. |12j . [2], p3], interprets 
the observed price as a noisy representation of an unobserved one, typically assumed to 
be a continuous ltd semi-martingale: in this framework many important results exist on 
robust estimation of the realized volatility, but these models seem not tractable for dealing 
with high frequency trading problems and stochastic control techniques, mainly because 
the state variables are latent rather than observed. 

(ii) The micro-to-macroscopic approach (e.g. [I], [3], [6], [1], [9]) uses point processes, in 
particular Hawkes processes, to describe the piecewise constant observed price, that moves 
on a discrete grid. In contrast with the macro-to-microscopic approach, these models do 
not rely on the arguable existence assumption of a fair or fundamental price, and focus 
on observable quantities, which makes the statistical estimation usually simpler. Moreover, 
these models are able to reproduce several well-known stylized facts on high frequency data, 
see e.g. [Z], 0: 

• Micro structure noise: high-frequency returns are extremely anticorrelated, leading to 
a short-term mean reversion effect, which is mechanically explicable by the structure 
of the limit order book. This effect manifests through an increase of the realized 
volatility estimator {signature plot) when the observation frequency decreases from 
large to fine scales. 

• Volatility clustering: markets alternates, independently of the intraday seasonality, 
between phases of high and low activity. 

• At large scale, the price process displays a diffusion behaviour. 

In this paper, we aim to provide a tractable model of tick-by-tick asset price for liquid 
assets in a limit order book with a constant bid-ask spread, in view of application to optimal 
high frequency trading problem, studied in the companion paper [10]. We start from a 
model-free description of the piecewise constant mid-price, i.e. the half of the best bid and 
best ask price, characterized by a marked point process (T n , J n ) n , where the timestamps 
(T n ) represent the jump times of the asset price associated to a counting process (Nt), and 
the marks (J n ) are the price increments. We then use a Markov renewal process (MRP) 
for modeling the marked point process. Markov renewal theory [13] is largely studied in 
reliability for describing failure of systems and machines, and one purpose of this paper 
is to show how it can be applied for market microstructure. By considering a suitable 
Markov chain modeling for (J n ), we are able to reproduce the mean-reversion of price 
returns, while allowing arbitrary jump size, i.e. the price can jump of more than one tick. 
Furthermore, the counting process (Nt), which may depend on the Markov chain, models 
the volatility clustering, i.e. the presence of spikes in the volatility of the stock price, and 
we also preserve a Brownian motion behaviour of the price process at macroscopic scales. 
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Our MRP model is a rather simple but robust model, easy to understand, simulate and 
estimate, both parametrically and non-parametrically, based on i.i.d. sample data. An 
important feature of the MRP approach is the semi Markov property, meaning that the 
price process can be embedded into a Markov system with few additional and observable 
state variables. This will ensure the tractability of the model for applications to market 
making and statistical arbitrage. 

The outline of the paper is the following. In Section [2j we describe the MRP model 
for the asset mid-price, and provides statistical estimation procedures. We show the simu- 
lation of the model, and the semi-Markov property of the price process. We also discuss 
the comparison of our model with Hawkes processes used for modeling asset prices and 
microstructure noise. Section [3] studies the diffusive limit of the asset price at macroscopic 
scale. In Section [4], we derive analytical formula for the mean signature plot, and compare 
with real data. Finally, Section 5 is devoted to a conclusion and extensions for future 
research. 



2 Semi-Markov model 

We describe the tick-by-tick fluctuation of a univariate stock price by means of a marked 
point process (T n , J n )neN- The increasing sequence (T n ) n represents the jump (tick) times 
of the asset price, while the marks sequence {J n )n valued in the finite set 

E = {+1, -1, . . . , +m, -m} C Z \ {0} 

represents the price increments. Positive (resp. negative) mark means that price jumps 
upwards (resp. downwards). The continuous-time price process is a piecewise constant, 
pure jump process, given by 

N t 

Pt = Po + ^Jn, t>0, (2.1) 
n=l 

where (Nt) is the counting process associated to the tick times (T n ) n , i.e. 

n 

N t = inf{n:^T fc <t} 
k=l 

Here, we normalized the tick size (the minimum variation of the price) to 1, and the asset 
price P is considered e.g. as the last quotation for the mid-price, i.e. the mean between 



the best-bid and best-ask price. Let us mention that the continuous time dynamics (2.1) is 
a model-free description of a piecewise-constant price process in a market microstructure. 
We decouple the modeling on one hand of the clustering of trading activity via the point 
process (Nt), and on the other hand of the microstructure noise (mean-reversion of price 
return) via the random marks (J n ) n . 

2.1 Price return modeling 

We write the price return as 

Jn = JnXn, n > 1, (2.2) 



3 



where J n := sign(J n ) valued in {+1,-1} indicates whether the price jumps upwards or 
downwards, and := \J n \ is the absolute size of the price increment. We consider that 
dependence of the price returns occurs through their direction, and we assume that only 
the current price direction will impact the next jump direction. Moreover, we assume that 
the absolute size of the price returns are independent and also independent of the direction 
of the jumps. Formally, this means that 

• (Jn)n is an irreducible Markov chain with probability transition matrix 

( l+a+ l—a+ \ 
l+ 2 a- J (2-3) 

with a + ,a_ G [—1,1). 

• (£n)n is an i.i.d. sequence valued in {1, . . . , m}, independent of ( J n ), with distribution 
law: Pi = P[£„ = i] G (0,1), i = l,...,m. 

In this case, (J n ) is an irreducible Markov chain with probability transition matrix given 
by: 



Q 



( PlQ ■ ■ ■ PmQ 



\ PlQ ■ ■ ■ PmQ 



(2.4) 



We could model in general (J n )n as a Markov chain with transition matrix Q = (qij) 
involving 2m(2m — 1) parameters, while under the above assumption, the matrix Q in ( |2.4| ) 
involves only m + 1 parameters to be estimated. Actually, on real data, we often observe 
that the number of consecutive downward and consecutive upward jumps are roughly equal. 
We shall then consider the symmetric case where 

a + = a_ =: a. (2-5) 

In this case, we have a nice interpretation of the parameter a G [—1, 1). 

Lemma 2.1 In the symmetric case, the invariant distribution of the Markov chain (*/n)n 

is tt = (212)) an ^ the invariant distribution of (J n ) n is it = (pi7r, . . . ,p m Tt). Moreoever, 
we have: 

a = corr„-(j n , J„_i), Vn > 1, (2.6) 

where coir n denotes the correlation under the stationary probability ¥ n starting from the 
initial distribution n. 

Proof. We easily check that under the symmetric case, nQ = n for tt = (1/2, 1/2), which 
means that fr is the invariant distribution of the Markov chain ( J n ) n . Consequently, it = 
(pifr, . . . ,p m ft) satisfies ttQ = tt, i.e. tt is the invariant distribution of (J n )n, and so under 
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IPVj (Jn)n (resp. (J n )n) is distributed according to 7r (resp. it). Therefore, E„-[j n ] = and 
Var^f/n] = 1. We also have for all n > 1, by definition of Q: 



IE ff [j n «/n-l] — E„ 



fl ^a-l) 2 -^(4-l) 2 



which proves the relation (2.6) for a. □ 

Lemma |2 . 1 1 provides a direct interpretation of the parameter a as the correlation between 
two consecutive price return directions. The case a = means that price returns are 
independent, while a < (resp. a > 0) corresponds to a mean-reversion (resp. trend) of 
price return. 

We also have another equivalent formulation of the Markov chain, whose proof is trivial 
and left to the reader. 

Lemma 2.2 In the symmetric case, the Markov chain (J n ) n can be written as: 

where (B n ) n is a sequence ofi.i.d. random variables with Bernoulli distribution on {+1, — 1}, 
and parameter (1 + a)/ '2, i.e. of mean¥,[B n ] = a. The price increment Markov chain [Jn)n 
can also be written in an explicit induction form as: 

Jn — Jn—lCm (^■'^) 

where (Cn)n is a sequence of i.i.d. random variables valued in E = {+1, —1, . . . , +m, — m}, 
and with distribution: P[Cn = k] = Pfc(l + sign(k)a)/2. 



The above Lemma is useful for an efficient estimation of a. Actually, by the strong law 
of large numbers, we have a consistent estimator of a: 

The variance of this estimator is known, equal to 1/n, so that this estimator is efficient and 
we have a confidence interval from the central limit theorem: 

V^(a (n) -a) ^> AA(0,1), as n -> oo. 

The estimated parameter for the chosen dataset is a = —87.5%, which shows as expected 
a strong anticorrelation of price returns. In the case of several tick jumps m > 1, the 
probability pi = P[£ n = i] may be estimated from the classical empirical frequency: 

p i ^L 1 ^' » = l,...,m 

fc=i 
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2.2 Tick times modeling 

In order to describe volatility clustering, we look for a counting process (Nt) with an 
intensity increasing every time the price jumps, and decaying with time. We propose a 
modeling via Markov renewal process. 

Let us denote by S n = T n — T n -\, n > 1, the inter- arrival times associated to (Nt). We 
assume that conditionally on the jump marks (J n )n, (S n )n is an independent sequence of 
positive random times, with distribution depending on the current and next jump mark: 

Fij(t) = F[S n+1 <t\J n = i, J n+1 = j], (i,j) e E. 

We then say that (T n , J n ) n is a Markov Renewal Process (MRP) with transition kernel: 

FVn+i = j, S n+ i <t\J n = i] = qijFij(t), (i,j) £ E. 

In the particular case where F^ does not depend on the point process (Nt) is indepen- 
dent of (J n ), and called a renewal process, Moreover, if F is the exponential distribution, 
N is a Poisson process. Here, we allow in general dependency between jump marks and 
renewal times, and we refer to the symmetric case when Fij depends only on the sign of ij, 
by setting: 

F + (t) = Fij(t), if ij>0, F_(t) = Fij(t), if ij<0. (2.8) 

In other words, F+ (resp. F-) is the distribution function of inter-arrival times given 
two consecutive jumps in the same (resp. opposite) direction, called trend (resp. mean- 
reverting) case. Let us also introduce the marked hazard function: 

h tj (t) = lim -F[t < S n+1 <t + S, J n +i =j\S n >t,J n = i], t>0, (2.9) 
54,0 

for i, j G E, which represents the instantaneous probability that there will be a jump with 
mark j, given that there were no jump during the elapsed time t, and the current mark is 
i. By assuming that the distributions F^ of the renewal times S n admit a density fij, we 
may write h{j as: 

/*'(*) 

hij(t) = Qij 1 _ ff. U\ =: 



where 



Hi(t) = F[S n+1 <t\J n = i] = ^qijF^t), 



is the conditional distribution of the renewal time in state i. In the symmetric case (2.5), 
(2.8), we have 

' 1 + sign(ij)cr 



, m / 1 + sign^ij Ja\ 
tiij(t) = pjy JA s i gn (y)(t), 
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with jump intensity 



where f± is the density of F±, and 

, . 1 + a „ . . 1 — a _ . . 

F(t) = —F+{t) + —F-(t). 

Markov renewal processes are used in many applications, especially in reliability. Cla- 
ssical examples of renewal distribution functions are the ones corresponding to the Gamma 
and Weibull distribution, with density given by: 

t /9-l e -t/0 ^ ^ fi(t\P-l 

where j3 > 0, and 6 > are the shape and scale parameters, and T (resp. r t ) is the Gamma 
(resp. lower incomplete Gamma) function: 

T(/3) = / s^e-'ds, T t {/3) = / sP^e-'ds, (2.11) 
Jo Jo 



fGam(t) = •— ^ f Wei (t) = tilLY ~ e -M 



2.3 Statistical procedures 

We design simple statistical procedures for the estimation of the distribution and jump 
intensities of the renewal times. Over an observation period, pick a subsample of i.i.d. 
data: 

{S k = T k - T k _i : k s.t. J fc _! =i, J k = j}, 

and set: 

hj = #{k s.t. J fc _i =i, J k = j} 
h = #{k s.t. J fc _i = i}, 

with cardinality respectively riij and rij. In the symmetric case, we also denote by 

I± = #{k s.t. sign(J k _iJ k ) = ±}, 

with cardinality n±. We describe both parametric and nonparametric estimation. 

• Parametric estimation 

We discuss the parametric estimation of the distribution F^ of the renewal times when 
considering Gamma or Weibull distributions with shape and scale parameters f3ij, 9ij. We 
can indeed consider the Maximum Likelihood Estimator (MLE) d%j), which are solution 
to the equations: 



m/ 3 I!M = i n (±Ys k )-±-Y\nS k 

1 lJ k=i k=i 



1 1 

A-S ni ., with S nij := —J2 Sk 
Pij 13 k=l 
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i j 


shape scale 


+1 +1 
-1 -1 
+1 -1 
-1 +1 


0.27651097 2187 
0.2806104 2565.371 
0.07442401 1606.308 
0.06840708 1508.155 



Table 1: Parameter estimates for the renewal times of a Gamma distribution 



There is no closed- form solution for which can be obtained numerically e.g. by Newton 
method. Alternatively, since the first two moments of the Gamma distribution S ~» T(ft, 9) 
are explicitly given in terms of the shape and scale parameters, namely: 

p _mf i__nsi (212) 

P ~ VarfS] ' 9 ~ Var[5]' 1 ' 



we can estimate ftij and 9%j by moment matching method, i.e. by replacing in (2.12) the 
mean and variance by their empirical estimators, which leads to: 



ft. 



n ijSn i:j 1 _ n-ijSriij 



)2 

k=l k=l 



We performed this parametric estimation method for the Euribor on the year 2010, 
from lOh to 14h, with one tick jump, and obtain the following estimates in Table [TJ 

We observed that the shape and scale parameters depend on the product ij rather than 
on i and j separately. In other words, the distribution of the renewal times are symmetric 
in the sense of ( |2.8| ). Hence, we performed again in Table [2] the parametric estimation 
(ft + , 9 + ) and (/?_,#_) for the shape and scale parameters by distinguishing only samples for 
ij = 1 (the trend case) and ij = —1 (the mean- reverting case). We also provide in Figure 
[l] graphical tests of the goodness-of-fit for our estimation results. The estimated value ft + , 
and ft- < 1 for the shape parameters, can be interpreted when considering the hazard rate 
function of S n +i given current and next marks: 

A±(t) := lim]:Flt<S n+1 <t + 8\S n >t,sign{J n J n+l ) = ±} 

t >0, 



i-F ± (ty 

when F± admits a density f± (Notice that X± differs from X±). \+(t) (resp. X-(t)) is the 
instantaneous probability of price jump given that there were no jump during an elapsed 
time t, and the current and next jump are in the same (resp. opposite) direction. For the 
Gamma distribution with shape and scale parameters (ft, 9), the hazard rate function is 
given by: 

-l 

/ / \ 

1 



^Gam 



it) 



e 



T(ft)-r t/0 (fty 
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i j 


shape scale 


+ 1 (trend) 
-1 (mean- reverting) 


0.276225 2397.219 
0.07132677 1561.593 



Table 2: Parameter estimates in the trend and mean-reverting case 



and is decreasing in time if and only if the shape parameter (3 < 1, which means that 
the more the time passes, the less is the probability that an event occurs. Therefore, the 
estimated values /3+, and /3_ < 1 are consistent with the volatility clustering: when a jump 
occurs, the probability of another jump in a short period is high, but if the event does not 
take place soon, then the price is likely to stabilize. Using this modeling, a long period of 
constant price corresponds to a renewal time in distribution tail. On the contrary, since 
renewal times are likely to be small when (3 < 1, most of the jumps will be extremely close. 
We also notice that the parameters in the trend and mean-reverting case differ significantly. 
This can be explained as follows: trends, i.e. two consecutive jumps in the same direction, 
are extremely rare (recall that a ~ —90%), since, in order to take place, market orders 
either have to clear two walls of liquidity or there must be a big number of cancellations. 
Since these events are caused by specific market dynamic, it is not surprising that their 
renewal law differ from the mean-reverting case. 




o 2000 4000 mi oooo 10000 o 2000 mm som aooo 10000 12000 

Gamma 10.0713 0.00064 ) 0anrn(o#6 , 0.W0417 ) 



c 




2000 4000 6000 0000 10000 2000 4000 6000 B000 10000 12000 

Renewal times Renewal tiras 



Figure 1: QQ plot and histogram of F_(left) and F + (right). Euribor3m, 2010, 10h-14h. 
• Non parametric estimation 

Since the renewal times form an i.i.d. sequence, we are able to perform a non parametric 
estimation by kernel method of the density fij and the jump intensity. We recall the smooth 
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kernel method for the density estimation, and, in a similar way, we will give the one for hij 
and Xij. Let us start with the empirical histogram of the density, which is constructed as 
follows. For every collection of breaks 

{0 < t\ < ... < tu < oo}, Sr '■= tr+l ~ t r 

we bin the sample (Sk = Tk — Tk-i)k=i,...,n and define the empirical histogram of as: 

1 4tik a T. j I < .<?, <- , ,\ 



fir Tl>ii 



By the strong law of large numbers, when riij — > oo, this estimator converges to 

1 



[t r < Sk < £ 



r+lj 



which is a first order approximation of f%j(t r ). Anyway, this estimator depends on the 
choice of the binning, which has not necessarily small (nor equal) <5 r 's. The corresponding 
non parametric (smooth kernel) estimator of the density is given by a convolution method: 



ft® = W-Sk), 



keii 



where Ki(x) is a smoothing scaled kernel with bandwidth b. In our example we have chosen 
the Gaussian one, given by the density of the normal law of mean and variance b 2 . In 
practice, many softwares provide already optimized version of non-parametric estimation of 
the density, with automatic choice of the bandwidth and of the kernel (here Gaussian) . For 
example in R, this reduces to the function density, applied to the sample {Sk \ k G hj}- 
In the symmetric case, the histogram reduces to: 

rhst ,, , _ 1 £l±\t r <S k < tr+l} 

J± VM — x ' 

o r n± 

while the kernel estimator is given by: 

fl p (t) = — K b (t-Sk). 

Figure [2] shows the result obtained for the kernel estimation of f±(t), compared to the 
corresponding histogram. The non parametric estimation confirms the decreasing form of 
both of the densities, whose interpetation is that most of the jump of the stock price takes 
place in a short period of time (less than second), even though some renewal times can go 
to hours. 



We use similar technique to estimate the marked hazard function hij defined in (2.9) 
By writing 

UM ,. 1P[* < Sk <t + S, J k = j I Jk-i =»] 
hij(t) = hm 



5 P [Sk > s I Jk-i 
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Histogram and non-parametric estimation of f- 



Histogram and non-parametric estimation of f+ 



4000 6000 8000 1 0000 1 2000 
renewal times fs) 




renewal times (s) 



Figure 2: Nonparametric estimation of the densities /+ and /_ 



we can define the empirical histogram of hij as 



l e Iij | t r <S k < t r + 5} 



5 r #{k eIi\S k > t r } 
and the associated smooth kernel estimator by 

K?w = E K ^ - s *) 1 



#{keli | S k >t}' 
Notice that this nonparametric estimator is factorized as 

-Km 



m I \n. 



'1.1 



where -^f is the estimator of qij and ^h^it) is the kernel estimator of the jump intensity 



f ■ ■ (t) 

Mj(t) = i-H (t) • Thus, we can either estimate Ajj and multiply by the estimator of to 
obtain hij or vice versa, obtaining the same estimators. In the symmetric case, we have: 



h h f(t r ) 



while 



1 #{k € I±\t r < S k < tj. + 5} 
5 r #{k | S k > t r } 



kei± 



#{k\S k >t}' 



Figure [3] shows the result obtained for X±, compared to the corresponding histogram. 
The interpretation is the following: immediately after a jump the price is in an unstable 
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condition, which will probably leads it to jump again soon. If it does not happen, the price 
gains in stability with time, and the probability of a jump becomes smaller and smaller. 
Moreover, due to mean-reversion of price returns, the intensity of consecutive jumps in the 
opposite direction is larger than in the same direction, which explains the higher value of 
A- compared to A+. 



Non parametric estimation of lambda- Non parametric estimation of lambda* 




2.8e-05 0.35 1.4 2.9 5 7.7 12 18 26 39 60 93 170 440 0.07 4.9 11 22 35 60 88 130 190 310 480 730 1500 5800 

s s 



Figure 3: Nonparametric estimation of the jump intensities A + and A_ 



2.4 Price simulation 



The simulation of the price process (2.1) with a Markov renewal model is quite simple, 
much easier than any other point process modeling microstructure noise. This would allow 
the user, even when Nt is better fit by a more complex point process, to have a quick proxy 
for its price simulation. Choose a starting price Pq at initial time To. 

Initialization step: • Set Po = Pq 

• draw Jo from initial (e.g. stationary) law ir 

Inductive step: k — 1 — > k (next price and next timestamp) 

• Draw Jfc according to the probability transition, and set P k = Pk-i + Jk- 

• Draw S k ~» Fj k _ lJk , and set T k = T fc _i + S k - 

Once (Tfc,i\)fc g N is known, the price process is given by the piecewise constant process: 

Pt = Pk, T k <t<T k+1 . 

We show in Figure [4] some simulated trajectories of the price process based on the estimated 
parameters of the Euribor3m. 
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Figure 4: Simulation of 9 trajectories of the Euribor with estimated parameters on empirical data 



2.5 Semi-Markov property 

Let us define the pure jump process, 



(2.13) 



which represents the last price increment. Then, / is a semi-Markov process in the sense 
that the pair (It, St) is a Markov process, where 



n+l, 



is the time spent from the last jump price. Moreover, the price process (Pf) is embedded 
in a Markov process with three observable state variables: (Pt,It,St) is a Markov process 
with infinitesimal generator: 



£<p(p,i,s) 



dip 
ds 



which is written in the symmetric case (2.5) and (2.8) as 



dip 1 + a 



Os 



A- 00 



1 — a 



^2 p i ~ si S n (0i) -sign(i)i, 0) - cp(p, i, s)] , 



where X± is the jump intensity defined in (2.10). 
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2.6 Comparison with respect to Hawkes processes 

In a recent work [3] (see also [9]), the authors consider a tick- by-tick model for asset price 
by means of cross-exciting Hawkes processes. More precisely, the stock price jumps by one 
tick according to the dynamics: 

P t = N+-N t -, (2.14) 

where is a point process corresponding to the number of upward and downward jumps, 
with coupling stochastic intensities: 

Xf = Aoo+ f <p(t-u)dN?, (2.15) 

J — oo 

where A^ > is an exogenous constant intensity, and ip is a positive decay kernel. This 
model, which takes into account in the intensity the whole past of the price process, provides 
a better good-fit of intraday data than the MRP approach, and can be easily extended to 
the multivariate asset price case. However, it is limited to unitary jump size, and presents 
some drawbacks especially in view of applications to trading optimization problem. The 



price process in (2.14)-(2.15) is embedded into a Markovian framework only for the case of 



exponential decay kernel, i.e. for 

<P® = 7e _/3 % + (t), 

with < 7 < /3. In this case, the Markov price system consists of P together with the 
stochastic intensities (A + ,A _ ). But in contrast with the MRP approach, the additional 
state variables (A + ,A~) are not directly observable, and have to be computed from their 



dynamics (2.15), which requires a "precise" estimation of the three parameters Aqo, 7, and 
/3. In that sense, the MRP approach is more robust than the Hawkes approach when 
dealing with Markov optimization problem. Notice also that in the MRP approach, we can 
deal not only with parametric forms of the renewal distributions (which involves only two 
parameters for the usual Gamma and Weibull laws), but also with non parametric form 
of the renewal distributions, and so of the jump intensities. Simulation and estimation in 
MRP model are simple since they are essentially based on (conditional) i.i.d. sequence of 
random variables, with the counterpart that MRP model can not reproduce the correlation 
between inter-arrival jump times as in Hawkes model. 

We finally mention that one can use a combination of Hawkes process and semi Markov 
model by considering a counting process N t associated to the jump times (T n ), independent 
of the price return (J n ), and with stochastic intensity 

X t = Aso + 7 fe-W-^dNu, (2.16) 
Jo 

with Aoo > 0, and < 7 < /3. In this case, the pair (It,Xt) is a Markov process, and 
(Pt,It,Xt) is a Markov process with infinitesimal generator: 

C<p(p,i,£) = p(n-£)^ + l^q ij [<p(jp + j,j,t + 'y)-<p(p,i,£)], 
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which is written in the symmetric case (2.5) as: 

dip 



Cip(p,i,£) = P(n 



+ 



+ 



1 + a 



1 — a 



^2pj[(f(p + sign(i)j, sign(i)j, £ + 7) - tp(p, i, £)] 



J^PjMP ~ sign(i)j, -sign(i)j, £ + 7) - <p(p, 



3 Scaling limit 

We now study the large scale limit of the price process constructed from our Markov renewal 



model. We consider the symmetric case (|2.5|) and (|2.8|), and denote by 

F(t) 



1 j Ct , 1 1 — OL „ , . 



the distribution function of the sojourn time S n . We assume that the mean sojourn time 
is finite: 



00 1 
t dF(t) < 00, and we set A = — . 
M 



By classical regenerative arguments, it is known (see that the Markov renewal process 
obeys a strong law of large numbers, which means the long run stability of price process: 

Pt 

— — > c, a.s. 
t 

when t goes to infinity, with a limiting constant c given by: 



i,jeE 



where E = {1, —1, . . . ,m, —m}, puj = J °° tdFij(t), Q = (qij) is the transition matrix (2.4) 
of the embedded Markov chain (J n ) n , and tt = (7Tj) is the invariant distribution of {J n ) n . 



In the symmetric case (2.5) and (2.8), we have qij = P\i\0- + sign(j)a)/2, 7Tj = p\i\/2, and 
so c = 0. We next define the normalized price process: 



Pi 



(T) 



PtT 



t€[0, 1], 



and address the macroscopic limit of at large scale limit T — > 00. From the func- 
tional central limit theorem for Markov renewal process, we obtain the large scale diffusive 
behavior of price process. 



Proposition 3.1 



lim P( T ) 

T->-oo 



(d) 
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where W = (Wt)te[o,i] ^ s a standard Brownian motion, and a^, the macroscopic variance, 
is explicitly given by 



a 



A 



Var(£ n )+(E[£ n ]) 2 ^ +n 



a 



(3.1) 



Proof. From [TT], we know that a functional central limit theorem holds for Markov 
renewal process so that 



p(T) 



(d) 



when T goes to infinity (here —I means convergence in distribution), where a x is given by 
cr 2 = Act 2 with A = and 



a 2 



with 



H 



J + 9j ~ 9i 
g=(gi)i = (Jam-Q + n)-^ 
b = (pi), h = ^Qijj, 



and II is the matrix defined by = iTj. In the symmetric case (2.5), a straightforward 
calculation shows that 



and then 



Thus, 



sign(i)a^2jpj = sign(i)aE[£ n ], 
3=1 



a 



sign(«)- E[£ n ], i e E = {1,-1, .. .,m,-m}. 

1 — a 



j if ij > 

i + ^E^J ifj>0,i<0 
J-i^ntn] if i < 0,i > 0. 



Therefore, a direct calculation yields 



2a 



(7 



£A j + i ^(ek„]) 2 

Var(£ n ) + (E[£ n ]) 2 ^. 



□ 
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4 Mean Signature plot 



In this section, we aim to provide through our Markov renewal model a quantitative justifi- 
cation of the signature plot effect, described in the introduction. We consider the symmetric 
and stationary case, i.e.: 

(H) 



(i) The price return (J n ) n is given by (2.2) with a probability transition matrix Q for 
{J n )n = (sign( J n )) n in the form: 



Q 




for some a G [—1, 1) . 



(ii) (Nt) is a delayed renewal process: for n > 1, S n has distribution function F (indepen- 
dent of i, j), with finite mean p, = J °° t dF(t) =: 1/A < oo, and finite second moment, 
and Si has a distribution with density (1 — F(t))/p,. 

It is known that the process N is stationary under (H)(ii) (see [8]), and so the price 
process is also stationary under the stationary probability ¥ n , i.e. the distribution of 
Pt+ T — Pt does not depend on t but only on increment time r. In this case, the empirical 
mean signature plot is written as: 



V(r) 



E 7T [(P lT -P ( 



ir<T 



^E«[(P T - P ) 2 ] 



(i-l)r) 



(4.1) 



Notice that if Pt = crWt, where Wt is a Brownian motion, then V is a flat function: 
V(t) = a 2 , while it is well known that on real data V is a decreasing function on r with 
finite limit when r — > oo. This is mainly due to the anticorrelation of returns: on a short 
time-step the signature plot captures fluctuations due to returns that, on a longer time- 
steps, mutually cancel. We obtain the closed-form expression for the mean signature plot, 
and give some qualitative properties about the impact of price returns autocorrelation. The 
following results are proved in Appendix. 



Proposition 4.1 Under (H), we have: 



V(r) 



<? 2 oo + 



-2a(E[£ n ]) V " G a (r) 



1 — a / (1 — a)r 



where o is the macroscopic variance given in (3.1), and G a (t) := E[a *] is given via its 
Laplace- Stieltjes transform: 



G a (s) = l-A(l-a) 



s(l-aF(s))' 



a / 0, 



(4.2) 
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F(s) := f™e- st dF(t). Alternatively, G a is given directly by the integral form: 

G a {t) = i_A(-^)(t-(l-a)jf Q° a (u)du), (4.3) 

oo 

where Qa(t) = a n F*^ (t) , and F*^ is the n-fold convolution of the distribution func- 

n=0 

tionF, i.e. F< n \t) = Jl F< n ~ l \t - u)dF{u) , F*W> = 1. 

Corollary 4.1 Under (H), we obtain the asymptotic behavior of the mean signature plot: 

V(oo) : = lim V(t) = , 

T— ^OO 

V(Q+) := limy(r) = AE[$. 

riO+ 

Moreover, 

a(V(0 + ) -V(oo)) < 0. 

Remark 4.1 In the case of renewal process where F is the distribution function of the 
Gamma law with shape j3 and scale 6, it is known that F*^ is the distribution function a 
the Gamma law with shape n/3 and scale 9, and so: 

F *(n) (t) _ T t/e(n(3) 

where T is the Gamma function, and Tt is the lower incomplete Gamma functions defined in 



(2.11). Plugging into (4.3), we obtain an explicit integral expression of the mean signature 



plot, which is computed numerically by avoiding the inversion of the Laplace transform 



(4.2). Notice that in the special case of Poisson process for N, i.e. F is the exponential 



distribution of rate A, the function G a is explicitly given by : G a (t) = e x ^ 

Remark 4.2 The term equal to the limit of the mean signature plot when time step 
observation r goes to infinity, corresponds to the macroscopic variance, and t^(0 + ) = AE[£^] 
is the microstructural variance. Notice that while increases with the price returns auto- 
correlation a, the limiting term V(0 + ) does not depend on a, and the mean signature plot 
is flat if and only if price returns are independent, i.e. a = 0. In the case of mean-reversion 
(a < 0), V(0 + ) > cr^, while in the trend case (a > 0), we have: t^(0 + ) < cr^. We display 
in Figures [5] plot example of the mean signature plot function for a Gamma distribution 
when varying a. We also compare in Figure [6] the signature plot obtained from empirical 
data on the Euribor, the signature plot simulated in our model with estimated parameters, 
and the mean signature for a Gamma distribution with the estimated parameters. This 
example shows how the signature plot decreasing form is a consequence (rather coarse) 
of the microstructure noise, and that the shape parameters of the gamma law, which is 
responsible for the volatility cluster, is able to reproduce the convexity of the signature plot 
shape. 
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Gamma signature plot (shape = rate 0.5 ) Gamma signature plot (shape = rate 0.5 ) 




■ alpha: -0.75 

■ alpha: -0.5 

■ alpha: -0.25 



Figure 5: r — > V(t) when varying price return autocorrelation a. Left: a > 0. Right: a < 



Signature Plot 




500 1000 1500 TOO 2500 3000 

tau(s) 



Figure 6: Comparison of signature plot: empirical data, simulated and computed 

5 Conclusion and extensions 

In this paper we used a Markov renewal process (T n , J n ) n to describe the tick- by-tick 
evolution of the stock price and reproduce two important stylized facts, as the diffusive 
behavior and the decreasing shape of the mean signature plot towards the diffusive variance. 
Having in mind a direct application purpose, we decided to sacrifice the autocorrelation of 
inter-arrival times in order to have a fast and simple non-parametric estimation, perfect 
simulation and the suitable setup for a market making application, presented in a companion 
paper [10]. Aware of the model limits, and with an eye to statistical arbitrage, the next 
step is to extend the structure of the counting process Nt, for example to Hawkes process, 
to have a better fit, and the structure of the Markov chain (J n ) n to a longer memory binary 
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processes, able to recognize patterns. Another direction of study will be the extension to 
the multivariate asset price case. 



A Appendix: Mean signature plot 

The price process Pt is given by 

Pt 

where 



Pa + Jk = ^ L n lN t =n, 



(A.l) 



k=l 



n=0 



^2 J k, n> 1, L c 



0. 



fc=i 



Lemma A.l Under (H), we have 



<j- 



n- 



2a 1 - a n , .2 
(EKn]) , 



A \ — a 1 — a 



where a x is defined in (3.1). 

Proof. By writing that L n = L n _i + (, n J n , we have for n > 1 

E tt [-^n] = IEtt [L^_x + + 2£,nL n -lJ n ] 

n-1 

= IEtt [^n-i] + + 2(E[£ n ]) 2 ^E 7r [j fc J n ], 

jfc=l 
n-1 

= E^L^] + +2(E[e„]) 2 ^E 7r [j fc Jo] 

fc=i 



(A.2) 



where we used the fact that = 1, (£ n )n are i.i.d, and independent of (J n )n i n t ne second 
equality, and the stationarity of (J n )n hi the third equality. Now, from the Markov property 
of (Jk)k with probability transition matrix Q in (2.3) and (2.5), we have for any k > 1: 



IEtt [^l^fc-l] ^0 

'l + a\ f f 



1 -a 
2 



Jk-l'h 



Ett [Jjfc-iJo 



from which we obtain by induction: 



E w [j k J ] = a k . 



Plugging into (A.2), this gives 



2a 



By induction, we get the required relation for E n [L^A . 

Consequently, we obtain the following expression of the mean signature plot: 



□ 
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Proposition A.l Under (H), we have 



V(r) 



a 



2a 1 - G a (r) 
1 — a (1 — a)r 

where G a (t) := E>^] = EZo^A^t = n]. 



(E[£ n ]) 2 , r>0, 



(A.3) 



Proof. From (4.1) and (A.l), we see that the mean signature plot is written as 

: -^2E n [L 2 n ]F n [N T = n], 



V(r) 



n=0 



since the renewal process N is independent of the marks ( J n ). Together with the expression 
of E^fL^] in Lemma A.l, this yields 

a^E^iNr] 2a 1-G a {r) 



V(r) 



(mn]) 2 . 



A t 1 — a (1 — a)r 
Finally, since E„-[A^ T ] = At by stationarity of N, we get the required relation. 



□ 



We now focus on the finite variation function G a that we shall compute through its 
Laplace-Stieltjes transform: 



G a {s) := / e~ st dG a (t), s>0. 
Jo- 

We recall the convolution property for Laplace-Stieltjes transform 

G7H = G.H, 

where 

G * H(t) = f G(t-s)dH(s). 
Jo 

Let us consider the function Q a defined by: 

oo 

Q a (t) := ^a n P^ >n]. 



n=0 



Lemma A. 2 Under (H), we have for all a /0, 



G a 



l~-)Qa + 



a 



a 



A / l-F(s) 



Q a (s) = 1 + a- 

s\\-aF{s) 

Proof. 1. For any a ^ 0, t > 0, we have 

oo oo 

G a (t) = a n Wn[N t = n] = J2 «™0W > n] - P w [JV t > n + 1] 

n=0 n=0 
oo ^ oo 

= > n] - - a n F n [N t > n] - F % [N t > 0] 



(A.4) 

(A.5) 
(A.6) 



n=0 

Qa(t)--(Qa(t)-l 

a V 



n=0 
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which proves (A.5). 

2. Recall that for the delayed renewal process N, the first arrival time Si is distributed 
according to the distribution A with density A(l — F). Let us denote by iV the no-delayed 
renewal process, i.e. with all interarrival times S® = T® — T^_ 1 distributed according to F, 
and by the function defined similarly as in (A.4) with N replaced by N°. Then, 



Q a {t) = l + aJ2a n ^[N t >n + l] 

n=0 

oo 

= l + a^a"P w [5i + T° <t] 

n=0 

oo 

= l + aE^^a n P^[T° <t-Si\Si < t] 

n=Q 

= 1 + aE w [Q° a (t - Si)l Sl 



1 + aQl * A(t), 



(A.7) 



By taking the Laplace-Stieltjes transform in the above relation, and from the convolution 
property, we get 



Qa 



1 + aQ° a A. 



By same arguments as in (A.7) and (A. 8), we get Q° = 1 + aQ^F, and so 

oi 1 



1 - aF 



(A. 



(A.9) 



Now, from the relation A(t) = J* * A(l — F(u))du, and by taking Laplace-Stieltjes transform 
we get: 



A 

-( 

s 



A(s) = 'i{l-F(s)). (A.10) 
By substituting (A.9) and (A. 10) into ( A.8| ), we get the required relation (A. 6). □ 



From the relations (A.5)-(A.6) in the above Lemma, we immediately obtain the expres- 
sion (4.2) for the Laplace-Stieltjes transform G a . Let us now derive the alternative integral 
expression (4.3) for G a . 

Lemma A. 3 Under (H), we have for all a ^ 0: 



Q a (t) = l + Xt-X(l-a) Q° a (u)du 



(A.ll) 



with 



n=0 

and F*( n ^ is the n-fold convolution of the distribution function F . 
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Proof. We rewrite the expression (A. 6) Laplace-Stieltjes transform as 

T OO 

aX 



Q a {s) = 1 + — (l-F)^2(aFY 
.1 ^— > 



n=0 



aX 



n=0 



1 + ^J2a n ((F) n -(F) n+1 



(A.12) 



Let us now consider the function G° a (resp. Q°) defined similarly as for G a (resp. Q a ) 
with N replaced by N°, the no-delayed renewal process with all interarrival times S® = 
T° — T^_ 1 distributed according to F. Then, 

OO 

G° a (t) = $>*(P ff [fl?>n]-P ff [JV?>n+l]) 

n=0 

oo 

= ^a"(p.[T n °<t]-P.[T n ° +1 <tf 

n=0 

oo 

= ^ - F< n+ ^)(t). 

n=0 

Therefore, the Laplace-Stieltjes transform of is written also as 

oo 

G °a = E Q "(^ n -^ n+1 )' 



n=0 



By defining the function I®(t) := J G^(u)du, we then see from (A.12) that 

Q a (s) = l + aXf°, 

and thus 

Q a (t) = 1 + aX I G° a (u)du. 
Jo 



(A.13) 



Finally, by same arguments as in (A.5), we have 

r*o — 



a 



1 



a 



and plugging into (A.13), we get the required result. 



□ 



By using (A.5) and (A.ll), we then obtain the integral expression (4.3) of the function 
G a as in Proposition 4.1 Finally, we derive the asymptotic behavior of the mean signature 
plot. 



Proposition A. 2 Under (H) ; we get: 

V{oo) := lim V(t) = 



V(0 + ) := limV(r) = XE[£], 



(A.14) 
(A.15) 



and 



V(0 + ) > V(po) if and only if a < 0. 
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Proof. By observing that the function G a is stricly bounded in r by 1 for any a S [—1,1), 
we easily obtain from the expression (A. 3) the limit for V(t) when r goes to infinity. 

On the other hand, by substituting the integral formula ( |4.3| > of the function G a into the 
expression (A.3) of the mean signature plot, we have: 



V{r) - 
Next, by noting that 

we deduce that 



a 



2a 
1 — a 



1 



n 



1 



I -a) 



Q° a (s)ds 



lim 

r^0+ T J 



Q° a (s)ds = Q° a (0) 



lim V(t) 



a 



1 



E[f„]) , 



o 



(A.16) 



which gives ( A. 15 ) from the expression (3.1 ) of a . Finally, we immediately see from ( A.16 ) 
that V(0 + ) > V(oo) if and only if a < 0. □ 
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